#include "csr_matrix.h"
#include <cstdio>
#include <algorithm>
#include <numeric>
#include <tuple>
#include <vector>
#include <mkl.h>
#include <suitesparse/amd.h>
#include <suitesparse/cholmod.h>
#include "scope_guard.hpp"
#include <fstream>

template <typename T>
static void realloc_array(T*& arr, std::size_t size) {
    if (arr) {
        mkl_free(arr);
    }
    if (size) {
        arr = static_cast<T*>(mkl_malloc(sizeof(T[size]), 64));
    }
    else {
        arr = nullptr;
    }
}

const CsrMatrix CSR_MATRIX_DEFAULT;

void SetupCsrMatrix(CsrMatrix* mat, int size, int max_nnz) {
    mat->size = size;
    mat->max_nnz = max_nnz;
    realloc_array(mat->row_ptr, size + 1);
    realloc_array(mat->col_idx, max_nnz);
    realloc_array(mat->value, max_nnz);
}

void CopyCsrMatrix(CsrMatrix* dst, const CsrMatrix* src) {
    const int n = src->size;
    const int nnz = GetCsrNonzeros(src);
    SetupCsrMatrix(dst, n, nnz);
    std::copy_n(src->row_ptr, n + 1, dst->row_ptr);
    std::copy_n(src->col_idx, nnz, dst->col_idx);
    cblas_dcopy(nnz, src->value, 1, dst->value, 1);
}

void CopyCsrMatrixValues(CsrMatrix* dst, const CsrMatrix* src) {
    const int nnz = GetCsrNonzeros(src);
    cblas_dcopy(nnz, src->value, 1, dst->value, 1);
}

void DestroyCsrMatrix(CsrMatrix* mat) {
    SetupCsrMatrix(mat, -1, 0);
}

#define MIN_DIAG 1e-12

static bool match(int i, const std::vector<std::vector<std::pair<int, double>>>& x, int* p, int* vis) {
    for (auto [j, a] : x[i]) {
        if (std::abs(a) > MIN_DIAG && !vis[j]) {
            vis[j] = 1;
            if (p[j] == -1 || match(p[j], x, p, vis)) {
                p[j] = i;
                return true;
            }
        }
    }
    return false;
}

static void max_match(const std::vector<std::vector<std::pair<int, double>>>& x, int* p) {
    int n = static_cast<int>(x.size());
    std::vector<int> vis(n);
    std::vector<int> ip(p, p + n);
    for (int i = 0; i < n; ++i) {
        std::fill(vis.begin(), vis.end(), 0);
        if (ip[i] == -1 && !match(i, x, p, vis.data())) {
            std::puts("error while matching");
            std::exit(-1);
        }
    }
}

void ReadCsrMatrixMM2(CsrMatrix* mat, const char* filename) {
	std::ifstream inputFs(filename);
	int nrows, nnz;


	inputFs >> nrows >> nrows >> nnz;
	
	SetupCsrMatrix(mat, nrows, nnz + nrows);
	

}

void ReadCsrMatrixMM3(CsrMatrix* mat, const char* filename) {
     int i;
     int nnz,nrows;
     std::ifstream inputFs(filename);
     int irow,icol;
     double val = 0;

     inputFs >> nrows >> nrows >> nnz;
     mat->row_ptr = new int[nrows];
     mat->col_idx = new int[nrows];
     mat->value = new double[nnz];
     mat->row_ptr[0] = 0;

    while (inputFs >> irow >> icol >> val) {
        irow--;
        icol--;        
    }
     
     inputFs.close();
}

void ReadCsrMatrixMM1(CsrMatrix* mat, const char* filename) {
     int i;
     int nnz,nrows;
     FILE *fid = fopen(filename,"r");/*open file*/
     int irow,icol;
     double val = 0;
     if (!fid){
         fprintf(stderr, "Error opening file: %s\n", filename);
         exit(1);
     }
     fscanf(fid, "%d %d %d\n", &nrows, &nrows, &nnz);

     SetupCsrMatrix(mat, nrows, nnz );/*Allocate memory*/
     for (i = 0; i < nrows; i++) {
		mat->row_ptr[i] = 0;
	}

     for (i=0; i < (nnz); i++)
     {
        fscanf(fid, "%d %d %lg\n", &irow, &icol, &val);
              irow--; icol--;
              mat->col_idx[i] = icol;/*Column index*/
              mat->value[i] = val;/*value*/
              mat->row_ptr[irow + 1] = i + 1;/*Offset array*/
     }
     fclose(fid);
}


#undef MIN_DIAG

void CsrMatVec(const CsrMatrix* mat, const double* vec, double* sol) {
    std::fill_n(sol, mat->size, 0);
    for (int row = 0; row < mat->size; ++row) {
        for (int i = mat->row_ptr[row]; i < mat->row_ptr[row+1]; ++i) {
            int col = mat->col_idx[i];
            sol[row] += mat->value[i] * vec[col];
        }
    }
}

int CsrAmdOrder(const CsrMatrix* mat, int* p, int* ip) {
    int n = mat->size;
    //int nnz = GetCsrNonzeros(mat);
    //std::vector<int> ap(n + 2, 0);
    //std::vector<int> ai(nnz);
    //for (int i = 0; i < nnz; ++i) {
        //ap[mat->col_idx[i]+2]++;
    //}
    //std::partial_sum(ap.begin(), ap.end(), ap.begin());
    //for (int i = 0; i < n; ++i) {
        //for (int j = mat->row_ptr[i]; j < mat->row_ptr[i+1]; ++j) {
            //ai[ap[mat->col_idx[j]+1]++] = i;
        //}
    //}
    //auto ret = amd_order(mat->size, ap.data(), ai.data(), p, nullptr, nullptr);
    auto ret = amd_order(mat->size, mat->row_ptr, mat->col_idx, p, nullptr, nullptr);
    if (ip) {
        for (int i = 0; i < n; ++i) {
            ip[p[i]] = i;
        }
    }
    return -(ret == AMD_INVALID || ret == AMD_OUT_OF_MEMORY);
}

int CsrNdOrder(const CsrMatrix* mat, int* p, int* ip) {
    cholmod_common cm;
    cholmod_start(&cm);
    int n = mat->size;
    int nnz = GetCsrNonzeros(mat);
    auto sp = cholmod_spzeros(n, n, nnz, CHOLMOD_PATTERN, &cm);
    std::copy_n(mat->row_ptr, n + 1, static_cast<int*>(sp->p));
    std::copy_n(mat->col_idx, nnz, static_cast<int*>(sp->i));

    //auto spt = cholmod_transpose(sp, 0, &cm);
    //auto spa = cholmod_add(sp, spt, nullptr, nullptr, false, true, &cm);
    //cholmod_free_sparse(&sp, &cm);
    //cholmod_free_sparse(&spt, &cm);
    //sp = spa;

    sp->stype = 1;
    std::vector<int> tmp(n + n);
    cholmod_nested_dissection(sp, nullptr, 0, p, tmp.data(), tmp.data() + n, &cm);
    //cholmod_metis(sp, nullptr, 0, false, p, &cm);
    auto ret = (cm.status == CHOLMOD_OK);
    cholmod_free_sparse(&sp, &cm);
    cholmod_finish(&cm);
    if (ip) {
        for (int i = 0; i < n; ++i) {
            ip[p[i]] = i;
        }
    }
    return ret;
}
